home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / isobuild / normal.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-02-11  |  16.2 KB  |  482 lines
  1.  
  2. /* normal.c              -Brian Tierney, LBL        1/91  */
  3.  
  4. /*  for use with isobuild:  computes normal vectors
  5.  *
  6.  *  There are 2 methods: 1 is to compute all normals ahead of time
  7.  *  and store them in an array.  This method is MUCH faster for the
  8.  *  server version, but uses LOTS of memory.  The other method is
  9.  *  to compute them as needed.
  10.  *
  11.  *  2 types of normal vectors are supported, and which type is controlled
  12.  *  by some #defines in isobuild.h
  13.  *
  14.  *  These types are:
  15.  *     1) default: normals will be in the range -1.0 to 1.0
  16.  *
  17.  *     2) #define SHORT_NORMALS :  normals will be in the range
  18.  *         -32767 to 32767 (short ints)  (or equal to the
  19.  *            range of the input data if input is type float)
  20.  *
  21.  *    It is not allowed to define both of these at the same time!
  22.  */
  23.  
  24. /* $Id: normal.c,v 1.4 1992/01/31 02:05:45 tierney Exp $ */
  25.  
  26. /* $Log: normal.c,v $
  27.  * Revision 1.4  1992/01/31  02:05:45  tierney
  28.  * *** empty log message ***
  29.  *
  30.  * Revision 1.3  1992/01/30  20:05:03  davidr
  31.  * prior to Brian's changes
  32.  *
  33.  * Revision 1.2  1991/12/19  01:42:20  davidr
  34.  * added RCS identification markers
  35.  * */
  36.  
  37. static char rcsid[] = "$Id: normal.c,v 1.4 1992/01/31 02:05:45 tierney Exp $" ;
  38.  
  39. #include "isobuild.h"
  40.  
  41.  
  42. /********************************************************************/
  43. void
  44. compute_normal_array()
  45. {                /* pre-compute all normals */
  46.     register int i, j, k;
  47.     NORMAL_VECT calc_normal();
  48.  
  49.     Status("Pre-Computing normals..");
  50.  
  51.     memset(gradient_curr_slice[0], 0, sizeof(FLOAT_VECT) * xdim * ydim);
  52.     memset(gradient_next_slice[0], 0, sizeof(FLOAT_VECT) * xdim * ydim);
  53.  
  54.     for (i = 0; i < zmax; i++) {
  55.     if (i % 10 == 0)
  56.         fprintf(stderr, "Computing normals: slice %d...\n", i);
  57.  
  58.     for (j = 0; j < ymax; j++)
  59.         for (k = 0; k < xmax; k++)
  60.         normals[i][j][k] = calc_normal(k, j, i);
  61.     }
  62.  
  63.     Status("Normals computed.");
  64. }
  65.  
  66. /**************************** calc_normal ****************************/
  67.  
  68. NORMAL_VECT
  69. calc_normal(x, y, z)
  70.     int       x, y, z;        /* location in original data set */
  71. {
  72.     static int prev_z = 0;
  73.     FLOAT_VECT **tmp_ptr, fvect, get_gradient();
  74.     NORMAL_VECT norm;        /* normals (returned by this routine) */
  75.     NORMAL_VECT normalize_vect();
  76.     FLOAT_VECT compute_cube_center_normal();
  77.     NORMAL_VECT normal_roundoff();
  78.  
  79.  
  80.     if (DIVIDING && DO_CUBE_CENTER_NORMALS) {
  81.     if (z != prev_z) {
  82.         tmp_ptr = gradient_curr_slice;
  83.         gradient_curr_slice = gradient_next_slice;
  84.         gradient_next_slice = tmp_ptr;
  85.         memset(gradient_next_slice[0], 0, sizeof(FLOAT_VECT) * xdim * ydim);
  86.         prev_z = z;
  87.     }
  88.     /* only compute gradient if it hasn't already been computed */
  89.     /*
  90.      * Actually: this method will re-compute the gradient if the previous
  91.      * gradient value was zero, but this doesnt happen often
  92.      */
  93.     if (gradient_curr_slice[y][x].z == 0.)
  94.         gradient_curr_slice[y][x] = get_gradient(x, y, z);
  95.     if (gradient_curr_slice[y][x + 1].z == 0.)
  96.         gradient_curr_slice[y][x + 1] = get_gradient(x + 1, y, z);
  97.     if (gradient_curr_slice[y + 1][x].z == 0.)
  98.         gradient_curr_slice[y + 1][x] = get_gradient(x, y + 1, z);
  99.     if (gradient_curr_slice[y + 1][x + 1].z == 0.)
  100.         gradient_curr_slice[y + 1][x + 1] = get_gradient(x + 1, y + 1, z);
  101.     /* compute values for the next slice */
  102.     if (gradient_next_slice[y][x].z == 0.)
  103.         gradient_next_slice[y][x] = get_gradient(x, y, z + 1);
  104.     if (gradient_next_slice[y][x + 1].z == 0.)
  105.         gradient_next_slice[y][x + 1] = get_gradient(x + 1, y, z + 1);
  106.     /* these 2 will have never been computed already */
  107.     gradient_next_slice[y + 1][x] = get_gradient(x, y + 1, z + 1);
  108.     gradient_next_slice[y + 1][x + 1] = get_gradient(x + 1, y + 1, z + 1);
  109.  
  110.     fvect = compute_cube_center_normal(x, y, z);
  111.  
  112.     } else {
  113.     fvect = get_gradient(x, y, z);
  114.     }
  115.  
  116.     norm = normalize_vect(fvect);
  117.  
  118.     return (norm);
  119. }
  120.  
  121. /******************************************************/
  122. NORMAL_VECT
  123. normalize_vect(vect)
  124.     FLOAT_VECT vect;
  125. {
  126.     double    sum, mag;
  127.     NORMAL_VECT norm;
  128.  
  129.     sum = (vect.x * vect.x) + (vect.y * vect.y) + (vect.z * vect.z);
  130.     mag = sqrt(sum + 0.0000001);
  131.  
  132. #ifdef SHORT_NORMALS
  133.     norm.x = (short) ((vect.x / mag) * 32767.0);
  134.     norm.y = (short) ((vect.y / mag) * 32767.0);
  135.     norm.z = (short) ((vect.z / mag) * 32767.0);
  136. #else
  137.     norm.x = (Norm_type) (vect.x / mag);
  138.     norm.y = (Norm_type) (vect.y / mag);
  139.     norm.z = (Norm_type) (vect.z / mag);
  140. #endif
  141.  
  142.     return (norm);
  143. }
  144.  
  145. /******************************************************************/
  146. FLOAT_VECT
  147. get_gradient(x, y, z)
  148.     int       x, y, z;
  149. /* This routine computes the gradient at location x,y,z */
  150. {
  151.     /* compute gradients using center differences */
  152.  
  153.     FLOAT_VECT grad;
  154.     Data_type d1, d2, val;
  155.  
  156. /*#define THIN_OBJ_CHECK doesn't really seem to help */
  157.  
  158. #define AVERAGE(d1,d2) (((float)d1 - (float) d2 ) * .5)
  159. #define MINABS(A,B)  ((fabs(A)) < (fabs(B)) ? (A) : (B))
  160.  
  161. #ifdef FLOAT_INPUT
  162.     if (x == 0)            /* use edge value instead */
  163.     grad.x = (data[z][y][x + 1] - data[z][y][x]) / 2.;
  164.     else if (x >= xmax)
  165.     grad.x = (data[z][y][x] - data[z][y][x - 1]) / 2.;
  166.     else
  167.     grad.x = (data[z][y][x + 1] - data[z][y][x - 1]) / 2.;
  168.  
  169.     if (y == 0)
  170.     grad.y = (data[z][y + 1][x] - data[z][y][x]) / 2.;
  171.     else if (y >= ymax)
  172.     grad.y = (data[z][y][x] - data[z][y - 1][x]) / 2.;
  173.     else
  174.     grad.y = (data[z][y + 1][x] - data[z][y - 1][x]) / 2.;
  175.  
  176.     if (z == 0)
  177.     tmpz = (data[z + 1][y][x] - data[z][y][x]) / 2.;
  178.     else if (z >= zmax)
  179.     tmpz = (data[z][y][x] - data[z - 1][y][x]) / 2.;
  180.     else
  181.     tmpz = (data[z + 1][y][x] - data[z - 1][y][x]) / 2.;
  182.  
  183. #else
  184.  
  185.     /* cast data to float to maintain more precision */
  186.     val = data[z][y][x];
  187.  
  188.     if (x == 0) {
  189.     d1 = data[z][y][x + 1];
  190.     d2 = val;
  191.     } else if (x >= xmax) {
  192.     d1 = val;
  193.     d2 = data[z][y][x - 1];
  194.     } else {
  195.     d1 = data[z][y][x + 1];
  196.     d2 = data[z][y][x - 1];
  197.     }
  198. #ifdef THIN_OBJ_CHECK
  199.     if ((d1 > val && d2 > val) || (d1 < val && d2 < val)) {
  200.     /* checks for 'thin' objects and computes gradient differently */
  201.     grad.x = MINABS(AVERAGE(d1, val), AVERAGE(d2, val));
  202.     } else
  203. #endif
  204.     grad.x = AVERAGE(d1, d2);
  205.  
  206.     if (y == 0) {
  207.     d1 = data[z][y + 1][x];
  208.     d2 = val;
  209.     } else if (y >= ymax) {
  210.     d1 = val;
  211.     d2 = data[z][y - 1][x];
  212.     } else {
  213.     d1 = data[z][y + 1][x];
  214.     d2 = data[z][y - 1][x];
  215.     }
  216. #ifdef THIN_OBJ_CHECK
  217.     if ((d1 > val && d2 > val) || (d1 < val && d2 < val)) {
  218.     grad.y = MINABS(AVERAGE(d1, val), AVERAGE(d2, val));
  219.     } else
  220. #endif
  221.     grad.y = AVERAGE(d1, d2);
  222.  
  223.     if (z == 0) {
  224.     d1 = data[z + 1][y][x];
  225.     d2 = val;
  226.     } else if (z >= zmax) {
  227.     d1 = val;
  228.     d2 = data[z - 1][y][x];
  229.     } else {
  230.     d1 = data[z + 1][y][x];
  231.     d2 = data[z - 1][y][x];
  232.     }
  233. #ifdef THIN_OBJ_CHECK
  234.     if ((d1 > val && d2 > val) || (d1 < val && d2 < val)) {
  235.     grad.z = MINABS(AVERAGE(d1, val), AVERAGE(d2, val));
  236.     } else
  237. #endif
  238.     grad.z = AVERAGE(d1, d2);
  239.  
  240. #endif
  241.  
  242.     return (grad);
  243. }
  244.  
  245. /********************************************************************/
  246. FLOAT_VECT
  247. compute_cube_center_normal(x, y, z)
  248.     int       x, y, z;
  249. {
  250.     FLOAT_VECT norm_vect;
  251.  
  252.     if (z < zmax) {
  253.     norm_vect.x = (.125 * (gradient_curr_slice[y][x].x +
  254.                    gradient_curr_slice[y][x + 1].x +
  255.                    gradient_curr_slice[y + 1][x].x +
  256.                    gradient_curr_slice[y + 1][x + 1].x +
  257.                    gradient_next_slice[y][x].x +
  258.                    gradient_next_slice[y][x + 1].x +
  259.                    gradient_next_slice[y + 1][x].x +
  260.                    gradient_next_slice[y + 1][x + 1].x));
  261.  
  262.     norm_vect.y = (.125 * (gradient_curr_slice[y][x].y +
  263.                    gradient_curr_slice[y][x + 1].y +
  264.                    gradient_curr_slice[y + 1][x].y +
  265.                    gradient_curr_slice[y + 1][x + 1].y +
  266.                    gradient_next_slice[y][x].y +
  267.                    gradient_next_slice[y][x + 1].y +
  268.                    gradient_next_slice[y + 1][x].y +
  269.                    gradient_next_slice[y + 1][x + 1].y));
  270.  
  271.     norm_vect.z = (.125 * (gradient_curr_slice[y][x].z +
  272.                    gradient_curr_slice[y][x + 1].z +
  273.                    gradient_curr_slice[y + 1][x].z +
  274.                    gradient_curr_slice[y + 1][x + 1].z +
  275.                    gradient_next_slice[y][x].z +
  276.                    gradient_next_slice[y][x + 1].z +
  277.                    gradient_next_slice[y + 1][x].z +
  278.                    gradient_next_slice[y + 1][x + 1].z));
  279.     } else {
  280.     norm_vect.x = (.25 * (gradient_curr_slice[y][x].x +
  281.                   gradient_curr_slice[y][x + 1].x +
  282.                   gradient_curr_slice[y + 1][x].x +
  283.                   gradient_curr_slice[y + 1][x + 1].x));
  284.  
  285.     norm_vect.y = (.25 * (gradient_curr_slice[y][x].y +
  286.                   gradient_curr_slice[y][x + 1].y +
  287.                   gradient_curr_slice[y + 1][x].y +
  288.                   gradient_curr_slice[y + 1][x + 1].y));
  289.  
  290.     norm_vect.z = (.25 * (gradient_curr_slice[y][x].z +
  291.                   gradient_curr_slice[y][x + 1].z +
  292.                   gradient_curr_slice[y + 1][x].z +
  293.                   gradient_curr_slice[y + 1][x + 1].z));
  294.     }
  295.     return (norm_vect);
  296. }
  297.  
  298. /********************************************************************/
  299. FLOAT_VECT
  300. compute_cube_center_normal_old(x, y, z)
  301.     int       x, y, z;
  302. {
  303.     FLOAT_VECT norm_vect;
  304.  
  305.     if (z < zmax) {
  306.     norm_vect.x = (.125 * (normals[z][y][x].x + normals[z][y][x + 1].x +
  307.                normals[z][y + 1][x].x + normals[z][y + 1][x + 1].x +
  308.                normals[z + 1][y][x].x + normals[z + 1][y][x + 1].x +
  309.           normals[z + 1][y + 1][x].x + normals[z + 1][y + 1][x + 1].x));
  310.  
  311.     norm_vect.y = (.125 * (normals[z][y][x].y + normals[z][y][x + 1].y +
  312.                normals[z][y + 1][x].y + normals[z][y + 1][x + 1].y +
  313.                normals[z + 1][y][x].y + normals[z + 1][y][x + 1].y +
  314.           normals[z + 1][y + 1][x].y + normals[z + 1][y + 1][x + 1].y));
  315.  
  316.     norm_vect.z = (.125 * (normals[z][y][x].z + normals[z][y][x + 1].z +
  317.                normals[z][y + 1][x].z + normals[z][y + 1][x + 1].z +
  318.                normals[z + 1][y][x].z + normals[z + 1][y][x + 1].z +
  319.           normals[z + 1][y + 1][x].z + normals[z + 1][y + 1][x + 1].z));
  320.     } else {
  321.     norm_vect.x = (.25 * (normals[z][y][x].x + normals[z][y][x + 1].x +
  322.               normals[z][y + 1][x].x + normals[z][y + 1][x + 1].x));
  323.  
  324.     norm_vect.y = (.25 * (normals[z][y][x].y + normals[z][y][x + 1].y +
  325.               normals[z][y + 1][x].y + normals[z][y + 1][x + 1].y));
  326.  
  327.     norm_vect.z = (.25 * (normals[z][y][x].z + normals[z][y][x + 1].z +
  328.               normals[z][y + 1][x].z + normals[z][y + 1][x + 1].z));
  329.     }
  330.     return (norm_vect);
  331. }
  332.  
  333.  
  334. /***********************************************************/
  335.  
  336. norm_debug(x, y, z)
  337.     int       x, y, z;
  338. {
  339.     FLOAT_VECT grad1, grad2, grad3, grad4, grad5, grad6, grad7, grad8;
  340.     NORMAL_VECT norm_vect;
  341.     NORMAL_VECT normalize_vect();
  342.     FLOAT_VECT get_gradient();
  343.     int       sx, sy, sz;
  344.  
  345.     sx = x;
  346.     sy = y;
  347.     sz = z;
  348.  
  349.     /**************/
  350.     fprintf(stderr, "This point: (%d,%d,%d) = %d \n", x, y, z, data[z][y][x]);
  351.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  352.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  353.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  354.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  355.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  356.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  357.     grad1 = get_gradient(x, y, z);
  358.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  359.         grad1.x, grad1.y, grad1.z);
  360.  
  361.     /**************/
  362.     x = sx + 1;
  363.     y = sy;
  364.     fprintf(stderr, "next point, this slice: (%d,%d,%d) = %d \n", x, y, z,
  365.         data[z][y][x]);
  366.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  367.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  368.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  369.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  370.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  371.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  372.     grad2 = get_gradient(x, y, z);
  373.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  374.         grad2.x, grad2.y, grad2.z);
  375.  
  376.     /**************/
  377.     x = sx;
  378.     y = sy + 1;
  379.     fprintf(stderr, "next point, this slice: (%d,%d,%d) = %d \n", x, y, z,
  380.         data[z][y][x]);
  381.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  382.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  383.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  384.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  385.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  386.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  387.     grad3 = get_gradient(x, y, z);
  388.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  389.         grad3.x, grad3.y, grad3.z);
  390.  
  391.     /**************/
  392.     x = sx + 1;
  393.     y = sy + 1;
  394.     fprintf(stderr, "next point, this slice: (%d,%d,%d) = %d \n", x, y, z,
  395.         data[z][y][x]);
  396.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  397.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  398.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  399.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  400.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  401.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  402.     grad4 = get_gradient(x, y, z);
  403.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  404.         grad4.x, grad4.y, grad4.z);
  405.  
  406.  
  407.     /**************************************/
  408.     z = sz + 1;
  409.     x = sx;
  410.     y = sy;
  411.     fprintf(stderr, "This point, next slice: (%d,%d,%d) = %d \n",
  412.         x, y, z, data[z][y][x]);
  413.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  414.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  415.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  416.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  417.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  418.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  419.     grad5 = get_gradient(x, y, z);
  420.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  421.         grad5.x, grad5.y, grad5.z);
  422.  
  423.     /**************/
  424.     x = sx + 1;
  425.     y = sy;
  426.     fprintf(stderr, "next point, next slice: (%d,%d,%d) = %d \n", x, y, z,
  427.         data[z][y][x]);
  428.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  429.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  430.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  431.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  432.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  433.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  434.     grad6 = get_gradient(x, y, z);
  435.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  436.         grad6.x, grad6.y, grad6.z);
  437.  
  438.     /**************/
  439.     x = sx;
  440.     y = sy + 1;
  441.     fprintf(stderr, "next point, next slice: (%d,%d,%d) = %d \n", x, y, z,
  442.         data[z][y][x]);
  443.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  444.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  445.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  446.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  447.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  448.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  449.     grad7 = get_gradient(x, y, z);
  450.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  451.         grad7.x, grad7.y, grad7.z);
  452.  
  453.     /**************/
  454.     x = sx + 1;
  455.     y = sy + 1;
  456.     fprintf(stderr, "next point, next slice: (%d,%d,%d) = %d \n", x, y, z,
  457.         data[z][y][x]);
  458.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  459.         x + 1, y, z, data[z][y][x + 1], x - 1, y, z, data[z][y][x - 1]);
  460.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  461.         x, y + 1, z, data[z][y + 1][x], x, y - 1, z, data[z][y - 1][x]);
  462.     fprintf(stderr, "     (%d,%d,%d) = %d ; (%d,%d,%d) = %d \n",
  463.         x, y, z + 1, data[z + 1][y][x], x, y, z - 1, data[z - 1][y][x]);
  464.     grad8 = get_gradient(x, y, z);
  465.     fprintf(stderr, "  gradient: gx = %f, gy = %f, gz = %f \n",
  466.         grad8.x, grad8.y, grad8.z);
  467.  
  468.     norm_vect.x = .125 * (grad1.x + grad2.x + grad3.x + grad4.x + grad5.x +
  469.               grad6.x + grad7.x + grad8.x);
  470.     norm_vect.y = .125 * (grad1.y + grad2.y + grad3.y + grad4.y + grad5.y +
  471.               grad6.y + grad7.y + grad8.y);
  472.     norm_vect.z = .125 * (grad1.z + grad2.z + grad3.z + grad4.z + grad5.z +
  473.               grad6.z + grad7.z + grad8.z);
  474.  
  475.     fprintf(stderr, "Normal for location (%d,%d,%d) = (%f,%f,%f) \n",
  476.         sx, sy, sz, norm_vect.x, norm_vect.y, norm_vect.z);
  477.  
  478.     norm_vect = normalize_vect(norm_vect);
  479.     fprintf(stderr, "Normalized: location (%d,%d,%d) = (%f,%f,%f) \n",
  480.         sx, sy, sz, norm_vect.x, norm_vect.y, norm_vect.z);
  481. }
  482.